CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     计算每个小球的Mie系数，
C     Output:
C        aMie bMie: Mie coefficients
C        rsr:
C        rsi:
C        rsx:
C        px:
C        nmax : the maximun order among all the spheres
C     Input:
C        x
C        ref
C        NADD
C        eps
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine Mie_init(nL, aMie, bMie, atr0, btr0, atr1, btr1, rsr,
     $     rsi, rsx, px, ek, drot, nmax, uvmax, k, r0, fint, NADD, eps)
      implicit double precision(a-h, o-z)
      include 'gmm01f.par'
      parameter (nmp=np*(np+2),nmp0=(np+1)*(np+4)/2)
      parameter (NXMAX=3000)
      parameter (ni0=np*(np+1)*(2*np+1)/3+np*np)
      parameter (ng0=np*(2*np**3+10*np**2+19*np+5)/6)
      parameter (nrc=4*np*(np+1)*(np+2)/3+np)
      parameter (nij=nLp*(nLp-1)/2)
C
      integer u,v,u0,nmax(nLp), uvmax(nLp), ind(nLp)
C
      double precision x(nLp), r0(6,nLp), w1(np), w2(np), w3(np),w4(np),
     +     rsr(np,nLp),rsi(np,nLp),rsx(np,nLp),px(np,nLp),
     +     rsr0(NXMAX),rsi0(NXMAX),rsx0(NXMAX),px0(NXMAX),
     +     drot(nrc,nij), k, besj(0:2*np+1),besy(0:2*np+1)
C
      complex*16 aMie(nLp, np), bMie(nLp, np), ref(nLp), an(np), bn(np)
     $     ,ek(np,nij), atr0(ni0, nij), btr0(ni0, nij), atr1(ni0, nij),
     $     btr1(ni0, nij)

      common/rot/    bcof(0:np+2),dc(-np:np,0:nmp)
      common/fnr/    fnr(0:2*(np+2))

      common/pitau/  pi(nmp0),tau(nmp0)
      common/tran/   atr

      common/ig0/    iga0(ni0)
      common/g0/     ga0(ng0)
      common/cofmnv0/cof0(ni0)
      common/crot/   cofsr(nmp)

C     constants
      pih    = dacos(0.d0)
      twopi  = 4.d0*pih
      pione  = 2.d0*pih
      ci     = dcmplx(0.d0, 1.d0)
      cin    = dcmplx(0.d0,-1.d0)

C      
C     adjust the radius of sphere according to the incident wave number
C
      do i=1,nL
         x(i)  = k*r0(4,i)
         ref(i)= dcmplx(r0(5,i),r0(6,i))
      enddo

      do j = 1,np
         do i = 1,nL
            aMie(i,j) = 0.d0
            bMie(i,j) = 0.d0
         enddo
      enddo

c
c  calculating Mie-scattering coefficients for each spheres
c  the ratio method of Wang and van der Hulst is used in calculating
c  Riccati-Bessel functions [see Wang and van der Hulst, Appl. Opt.
c  30, 106 (1991), Xu, Gustafson, Giovane, Blum, and Tehranian,
c  Phys. Rev. E 60, 2347 (1999)]
c
      nmax0 = 1
      do i = 1,nL

C     第一个球不检查重复性
         if(i.eq.1) goto  12

C     如果前后两个球一样，则可以节省计算：直接拷贝前一个球的数据
         if(x(i).eq.x(i-1).and.ref(i).eq.ref(i-1)) then
            nmax(i)  =  nmax(i-1)
            uvmax(i) = uvmax(i-1)
            goto 15
         endif

 12      write(6,'(a,i3,a,f7.2)') 'sphere #',i,
     +        '   individual size parameter: ',x(i)

C     nmax(i) is calculated here!
         call abMiexud(x(i),ref(i),np,NXMAX,nmax(i),an,bn,NADD,
     +                   rsr0,rsi0,rsx0,px0,w1,w2,w3,w4,eps)

         if(nmax(i).gt.np) then
            write(6,*) ' Parameter np too small, must be >',nmax(i)
            write(6,*) ' Please change np in gmm01f.par, recompile,'
            write(6,*) '   then try again'
            stop
         endif

         uvmax(i) = nmax(i)*(nmax(i)+2)
         write(6,'(a,1x,i4)')
     +        ' Actual single-sphere expansion truncation:',nmax(i)
         do j = 1,nmax(i)
            rsr(j,i) = rsr0(j)
            rsi(j,i) = rsi0(j)
            rsx(j,i) = rsx0(j)
            px(j,i)  = px0(j)

            temp1 = an(j)
            temp2 = bn(j)
            if(j.eq.1.or.j.eq.nmax(i))
     +           write(6,'(i10,4e15.7)') j,temp1,
     +           dimag(an(j)),temp2,dimag(bn(j))
         enddo

C     3.计算得到所有小球的Mie系数
 15      do j=1,nmax(i)
            aMie(i,j) = an(j)
            bMie(i,j) = bn(j)
            rsr(j,i)  = rsr0(j)
            rsi(j,i)  = rsi0(j)
            rsx(j,i)  = rsx0(j)
            px(j,i)   = px0(j)
         enddo

C     4. 返回最大的展开阶数
         if(nmax(i).gt.nmax0) nmax0 = nmax(i)
      enddo

C     
C     Mie系数计算完毕，开始计算其他系数
C     
C     
C     the formulation used here for the calculation of Gaunt coefficients
C     can be found in Bruning and Lo, IEEE Trans. Anten. Prop. Ap-19, 378
C     (1971) and Xu, J. Comput. Appl. Math. 85, 53 (1997), J. Comput. Phys.
C     139, 137 (1998)
C     
      n0 = nmax0 + 2
      fnr(0)=0.d0
      do n = 1,2*n0
         fnr(n) = dsqrt(dble(n))
      enddo
      bcof(0)=1.d0
      do n = 0,n0-1
         bcof(n+1) =fnr(n+n+2)*fnr(n+n+1)*bcof(n)/fnr(n+1)/fnr(n+1)
      enddo
      call cofsrd(nmax0)
      call cofd0(nmax0)
      call cofnv0(nmax0)
      call gau0(nmax0)
C     
      do i=1,nL
         do j=i+1,nL
C
            ij = (j-1)*(j-2)/2 + j - i

            x0 = r0(1,i) - r0(1,j)
            y0 = r0(2,i) - r0(2,j)
            z0 = r0(3,i) - r0(3,j)
            call carsphd(x0,y0,z0,d,xt,sphi,cphi)

            temp = (r0(4,i) + r0(4,j))/d
            if(temp.le.fint) goto 16

C     Constants 1:
            ephi  = dcmplx(cphi,sphi)
            nlarge= max(nmax(i),nmax(j))
            do m = 1,nlarge
               ek(m,ij) = ephi**m
            enddo
            xd = k*d
            nbes = 2*nlarge + 1
            call besseljd(nbes,xd,besj)
            call besselyd(nbes,xd,besy)
C
C  subroutine rotcoef.f calculates (-1)^{m+u)d_{mu}^n where d_{mu}^n is
C  the "reduced rotation matrix elements"
C  rotcoef.f is originally written by Mackowski (taken from scsmfo1b.for
C  developed by Mackowski, Fuller, and Mishchenko)
C
            call rotcoef(xt,nlarge)
            irc=0
            do n=1,nlarge
               n1=n*(n+1)
               do u=-n,n
                  do m=-n,n
                     imn = n1  + m
                     irc = irc + 1
C     Constants 2:
                     drot(irc,ij) = dc(u,imn)
                  enddo
               enddo
            enddo
            itrc = 0
            nsmall=min(nmax(i),nmax(j))
C     
C  the formulation used here for the calculation of vector translation
C  coefficients are from Cruzan, Q. Appl. Math. 20, 33 (1962) and
C  Xu, J. Comput. Phys. 139, 137 (1998)
C 
c MR: The translation coefficients atr0, btr0 get you from
c MR: basis {N3(j), M3(j)} to basis {N1(k), M1(k)}
c MR: The translation coefficients atr1, btr1 get you from
c MR: basis {N3(j), M3(j)} to basis {N3(k), M3(k)}
            do m=-nsmall,nsmall
               n1=max(1,iabs(m))
               do n=n1,nlarge
                  do v=n1,nlarge
                     itrc=itrc+1
C     Constants 3.
                     call cofxuds0(nmax0,m,n,v,besj,besy,
     +                    atr0(itrc,ij),btr0(itrc,ij),
     +                    atr1(itrc,ij),btr1(itrc,ij))
                  enddo
               enddo
            enddo
 16         continue
          enddo
       enddo
       
       return
       end
